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Abstract 

A model and analysis of the human sleep/wake system is presented. The 

model is derived using the known neuronal groups, and their various pro- 
jections, involved with sleep and wake. Inherent in the derivation is the 
existence of a slow time scale associated with homeostatic regulation, and 
a faster time scale associated with the dynamics within the sleep phase. A 
significant feature of the model is that it does not contain a periodic forcing 
term, common in other models, reflecting the fact that sleep/wake is not de- 
pendent upon a diurnal stimulus. Once derived, the model is analyzed using 
a linearized stability analysis. We then use experimental data from normal 
sleep-wake systems and orexin knockout systems to verify the physiological 
validity of the equations. 

Keywords: Sleep/ Wake Cycle, circadian regulation, limit cycle, Hopf 
Bifurcation 



1. Introduction 

Most living organisms exhibit daily biological rhythms. One of the most 
curious and difficult to understand is the sleep/wake cycle. Even though we 
are very familiar with sleep/wake, the reasons for sleep are not well under- 
stood. What is intriguing about the sleep- wake cycle is the complex interplay 
between the arousal and sleep states. Part of this is due to the cychc alter- 
ation of the rapid eye movement (REM) and non-REM (NREM) stages of 
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sleep. The distribution of REM and NREM varies dramatically across the 
animal kingdom, to the point that a clear pattern has not emerged. One 
consequence of this is that an animal model for human sleep has not been 
found ^42]. Another complication is that sleep/wake is sometimes thought 
to be synonymous with the 24 hour circadian clock. This is not correct. 
Although sleep/wake is coupled to the circadian clock, it is regulated by a 
neural system distinct from the circadian timing mechanism |3l]. Sleep/wake 
are also heavily influenced by two other factors. One is homeostatic, where 
the need for sleep increases with the length of the wake period. The second is 
allostatic, where the sleep/wake cycle is modified due to external behavioral 
events • The circadian process is connected to the sleep- wake cycle, while 
a homeostatic process is inherent to the sleep-wake cycle. To be specific, it 
is now known that the sleep-wake cycle is independent of the diurnal stim- 
ulus, and is instead driven by homeostatic forces [12], [T71 [1]. Therefore, the 
inherent periodicity of the sleep/wake cycle is not due to a periodic forcing 
mechanism associated with the diurnal cycle but is an inherent characteristic 
of the system. The objective of this study is to derive, and then analyze, a 
model for the sleep/wake cycle that accounts for the circadian and homeo- 
static influences of sleep/wake, but not necessarily the allostatic affects. This 
will include both REM and NREM sleep, as well as the waking state and 
observations about the associated brainwave activity of all states. 

1.1. Biological Background 

A basic assumption is that the brain has two states, sleep and wake, with 
sharp transitions between them. There are multiple cell groups in the brain 
contributing to these states, but the monoaminergic nuclei (MN) have been 
identified as primarily responsible for promoting the awake state [32] . Sleep is 
primarily promoted by the ventrolateral preoptic nucleus (VLPO), which is a 
group of cells associated with generating NREM and REM. A consequence of 
this distinct structural organization is that although the states are coupled, 
they have different regulatory mechanisms in the brain. 

1.1.1. REM /NREM 

The REM and NREM phases of sleep have a well defined pattern (refer 
to Figure [2]). Regardless of the conditions of sleep onset, a person always 
begins in an NREM period of the cycle. After approximately 80-85 minutes 
the brain switches to the REM stage, and remains in this state for approxi- 
mately 5-10 minutes before switching back to NREM. This 90 minute pattern 
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continues through the night with the REM-on intervals becoming progress- 
ing longer, and the NREM-on intervals shorter. The last REM-on period 
is about 25 minutes. There is no probability of waking naturally during 
an NREM interval, as the NREM component occurs during the periods of 
slow wave sleep. It is postulated that this is due to NREM being the most 
mentally restorative of the two processes, requiring little brain function and 
unconsciousness. Natural wake up occurs towards the end of the last NREM 
period. [39] 

There are three different ways of triggering sleep onset: exciting sleep 
active neurons, inhibiting wake active neurons or a combination of the two. 
The same idea applies to the triggering of wake onset. Consequently, to 
investigate the source of sleep/wake brain function, it is necessary to examine 
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Figure 2: The basic partitions of the states of sleep and wake for humans. 



the neurotransmitter activity of the system. 

1.1.2. Location Specific Neurotransmitter Activity 

The chemical activity of the brain involves the secretion of hormones 
and neurotransmitters from specific sections of the brain. The sleep-wake 
cycle is initiated and maintained by multiple distinct neuronal regions in 
the brainstem. These regions are located in the hypothalamus and are cou- 
pled to each other via electrical and neurochemical signaling. Figures [3] 
and |4] provide a schematic of location specific neurotransmitter signal flow 
in terms of excitatory and inhibitory mechanisms. Figure [T] shows the ap- 
proximate locations of sleep and wake promoting areas in the human brain. 
Figure [5] provides a simple depiction of the switching between wake and 
sleep states and well as the switch between REM and NREM states, and 
the set of neurotransmitters associated with each, as well as their excita- 
tory (+) or inhibitory (-) effects. The major neurotransmitter players are 
GABA (gamma-aminobutyric acid), orexin (hypocretin), histamine, acetyl- 
choline, noradrenaline (norepinephrine), serotonin and dopamine. Each of 
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these transmitters are located throughout the following neuronal groups: 
basal forebrain (BF), raphe nuclei (RN), laterodorsal tegmentum (LDT), 
pedunculopontine tegmentum (PPT), locus coeruleus (LC), lateral hypotha- 
lamus (LH), tuberomammillary nucleus (TMN), ventral periaqueductal gray 
region (vPAG), and ventrolateral preoptic nucleus (VLPO). Adenosine has 
been proposed to be the homeostatic accumulator of the need for sleep |32]- 
Acetylcholine effects muscle action, learning and memory. Norepinephrine 
effects one's control of alertness and wakefulness. Dopamine effects move- 
ment, focus and learning. Serotonin regulates mood, appetite, sleep, arousal 
and pain threshold. GABA contributes to sleep control and appetite control. 
Also note that the sections of the brain we will be focusing on are the cere- 
bral cortex, hypothalamus and thalamus, which play key roles in homeostatic 
maintenance. In this way they are the most relevant to the control of the 
human sleep/wake system. 

Figure [6] shows one full period of the sleep-wake cycle model in terms of 
adenosine and GABA concentration levels. The point at which the curves 
cross indicate wake or sleep initializations, depending on whether one of the 
variables is increasing or decreasing. Wake Initialization occurs when GABA 
is decreasing and AD is increasing. Sleep initialization occurs when AD is 
decreasing and GABA is increasing. 
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Figure 3: Excitatory projections of neurotransmitter flow in the sleep/wake generating 
network. 



Adenosine, which is a neuromodulator, plays a major role in the sleep- 
wake process. It is a byproduct of adenosine-triphosphate (ATP), the main 
"energy currency" of the brain, whenever neurotransmitters are fired, adeno- 
sine is produced. Adenosine is also continuously converted back to ATP via 
the metabolic processes that use it as an energy source [37j. The adenosine 
concentration that we are most interested in is that of the basal forebrain. 
While inhibiting the expression of all wake active neurotransmitters and dis- 
inhibiting the expression of those found in the VLPO [7] , adenosine levels in 
the basal forebrain increase during prolonged wake and fall during recovery 
sleep, providing an intrinsic periodic force associated with waking jSHl EH]- 
The linear relationship between between adenosine accumulation and sleep 
propensity provides support for using adenosine as representative of sleep 
propensity. GABA in the VLPO can be viewed as the other major player in 
the sleep-wake system. GABAa receptors in the brain mediate the tonic (or 
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Figure 4: Inhibitory projections of neurotransmitter flow in the sleep/wake generating 
network 
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Figure 5: The cycle of NREM Sleep, REM Sleep and Wakefulness with the neurotrans- 
mitters characteristic of their stages. 



slow) conductance |10]. This slower type of signaling plays a key role in cell 
excitability and sleep homeostasis [10] • Neuron populations in the VLPO 
are inactive during the vigilance state, and have a high firing rate and level 
during NREM, especially during delta wave (slow wave) sleep |10]. Being 
the major inhibitory neurotransmitter in the central nervous system, GABA 
is a key player in many physical systems. 

Referring to the Table 1, we can see the detailed activity associated with 
each neurotransmitter and brain area, allowing for the development of a sys- 
tem of variables. As a simple example, OXlh is the variable associated with 
the chemical concentration of orexin in the lateral hypothalamus. Since we 
only discuss orexin in one location, the subscript LH can be dropped. Orexin 
expression is promoted by acetylcholine in the LDT/PPT and noradrenaline 
in the LC. Orexin expression is inhibited by adenosine, serotonin, GABA 
and noradrenaline. Orexin promotes the release of noradrenaline, serotonin, 
histamine, GABA in the basal forebrain, acetylcholine, and dopamine, as 
well as promoting its own release (autoexcitation). These dynamics indicate 
the need for a linear combination of acetylcholine, noradrenaline, serotonin, 
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Figure 6: One full period of the sleep- wake cycle model. 



dopamine, histamine, GABA, adenosine, and orexin to express the rate of 
change of OX concentration. 



2. Materials and Methods 

In order to test the functionahty of the model, we use numerical methods 
to determine how the system responds to a variety of perturbations. In 
the calculations to follow, MATLAB 2009a was used to run all the programs 
written for the Human Sleep/Wake system described here. However, a fourth 
order Runge-Kutta method was used for solving the system of equations, 
instead of using one of the MATLAB ode solvers because of the unique 
perturbations used in the sleep experiments described. 
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Figure 7: Components of the human sleep/wake cycle as well as their associate neuro- 
transmitters 



3. Theory and Calculation 

3.1. Previous Work 

A good survey of the early sleep/wake models can be found in jSS]- We 
will only comment here on those most relevant to our development. 

In a recent paper by Phillips and Robinson [28], a model was proposed 
in order to improve quantitative understanding of the human sleep/wake cy- 
cle, which is similar to the intentions of this study. They use the following 
three neuronal populations to model the human ascending arousal system: 
monoaminergic (MA) nuclei, acetylcholine related (ACh) nuclei, and popu- 
lations in the VLPO. The VLPO is then split into circadian and homeostatic 
drives, which are represented by periodic forcing functions. Each population 
has its dynamics modeled by a sum of mean cell body potentials and mean 
neuronal firing rates subjected to the periodic forcing of the VLPO. While 
their model consistently predicts sleep/wake cycle dynamics observed in pre- 
vious studies, it is important to note that the homeostatic and circadian 
drives are represented as periodic forcing rather than dynamics inherent to 



10 



the physical system. 

[To] developed a model of the network of neurons controlling behavioral 
state transitions in mice, while also identifying special features of human 
and mouse sleep. It is similar to the Phillips and Robinson model (as well 
as the model being developed here) in that it defines three main processes 
being modeled: Sleep activity, Wake activity, and REM activity. Each pop- 
ulation is defined by a set of three equations, determined by a wide range in 
timescales. It incorporates homeostatic REM and NREM sleep drive and uses 
known sleep active neurotransmitters to determine parameters. It then uti- 
lizes extensive experimental data gathered from mouse sleep/wake states to 
verify the accuracy of the model. A limitation of their approach concerns the 
separation of vigilance states, which are characterized in a non-quantitative 
fashion. Similarly, "activity states" are used rather than specific neurotrans- 
mitter concentration and/or neuron firing data. 

Taking a somewhat different viewpoint, [31] use the "flip-flop" dynamics 
proposed in [32] to construct a biologically-based mathematical model that 
accounts for several features of the human sleep/wake cycle. This differs 
from our model in several ways, but two examples are that they assume, 
incorrectly, that the SCN projects directly to the VLPO, and an external 
periodic forcing mechanism is used to drive the model. 

In another recent paper by [9], a model of the rat sleep/wake regula- 
tory network is formed utilizing similar techniques as those used to model 
a mouse sleep/wake regulatory network in her aforementioned work [TO]. In 
this model, neurotransmitter concentrations are coupled to neuronal group 
firing rates to create a more complete view of the hypothalamic sleep-wake 
regulatory network. The model is composed of three parts: wake/sleep active 
neuron population firing rates, wake/sleep active neurotransmitter concentra- 
tions and a homeostatic regulator. All physical parameters were determined 
using extensive experimental data. This paper emphasizes the importance of 
utilizing neuron firing rates as well as neurochemical concentrations to form a 
physiologically accurate model of sleep/wake regulation. A downside of this 
model lies in the use of steady state saturation functions and Heaviside func- 
tions. Rather than using inherent neurotransmitter dynamics, the functions 
of neurotransmitter release are approximated using sigmoidal functions. The 
homeostatic drive is a physically unquantifiable force composed of Heaviside 
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switching mechanisms. 



3.2. Model Formulation 

In this model for the sleep/wake cycle, we are considering three principal 
states. They are the wake state iW), the NREM sleep state {S) and the 
REM sleep state {R). Table 1 provides a framework for the detailed interac- 
tions of all the neurotransmitters involved in the sleep-wake system. Given 
the number of variables involved, before deriving the equations it is worth 
describing the dynamics in more qualitative terms. 

The sleep-specific and wake specific neurotransmitters operate on a slow 
time scale, measured in terms of hours, while the REM and NREM specific 
neurotransmitters operate on a fast time scale system., measured in terms of 
minutes. The two neurochemicals associated with the slow time scale are 
adenosine in the BP and GABA in the VLPO. We use these two concentra- 
tions as variables indicating promotion of wake and sleep. High adenosine 
concentrations in the BP promote wakefulness, while high GABA concentra- 
tions and fast firing rates in the VLPO promote sleep. Adenosine's slow build 
up and decay is associated with wakefulness, while GABA in the VLPO is as- 
sociated with sleep. These observations are the basis for the hypothesis that 
adenosine is the principal component in determining wake, and GABAyipo 
is the principal component in promoting sleep [lOl [301 EH]. We also know 
from these references that AD disinhibits the VLPO during slow wave sleep, 
providing a direct connection between the building adenosine sleep pressure 
during waking and the gradual relief of this pressure during sleep. The other 
neurotransmitters are assumed to be associated with the fast time scale dy- 
namics. Adenosine and GABAvlpo fluctuate with a period of 24 hours, 
while the other neurotransmitters fluctuate with a period of several minutes. 

An important inherent quality of the sleep-wake cycle is that it should 
maintain it's state transitions in the absence of external periodic forcing. 
Mathematically, this indicates the need for the system to have a stable (at- 
tracting) limit cycle. This solution should also be stable to certain perturba- 
tions, such as those that might occur from genetic mutations causing the pro- 
duction of too much or too little of a neurotransmitter, or lacking the proper 
neuron groups to produce a neurotransmitter at all. Perturbations may also 
occur via the removal of specific neuron groups, blockage of neurotransmitter 
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bonding sites, and external forcing. It is well known from experimental data 
that perturbations to neurotransmitters don't alter the two-state nature of 
the sleep-wake system [SOlEaiaEIllIHlEaEIlEOlElIIIlEaElIH As 
an example, cholinergic agonists applied to the LC in cats leads to increased 
REM sleep, but the sleep-wake cycle as a whole remains periodic [S]. There- 
fore one of the central goals will be to locate the area in the phase plane of 
GABAvLPO and AD where a limit cycle exists and is stable. 

3.2.1. Kinetics 

For each neurotransmitter involved in the sleep-wake system, we iden- 
tify a location specific variable representing the effective neurotransmitter 
concentration. We assume first order kinetics for the neurotransmitters op- 
erating on the fast time scale. The fast time system may be formed from the 
observations compiled about the inhibitory and excitatory nature of the neu- 
rotransmitters in Table 1 combined with the dynamics arising in a synapse. 
If Xi is representative of a neurotransmitter in region i, we have 

dX- 

= Release(t) — Uptake(t) — Degradation(t). (3.1) 

Note that inhibitory effects are considered pre-synaptic and are therefore 
controlling flow from the neuron into the synaptic gap. Despite the fact that 
inhibitory and excitatory effects can occur pre- or post-synaptically, for sim- 
plicity we will assume the blocking or binding process of inhibition as only 
controlling flow into the gap. The Release function contains rate parame- 
ters and neurotransmitter variables that excite and inhibit the expression 
of Xi. The Uptake function depends on Xi and its associated vesicular up- 
take transporter, if existent. The Degradation function depends on Xi and 
an associated rate of decay. The uptake transporters are modeled using 
known logistic dynamics from neurotransmitter uptake transporter micro- 
dialysis measurements [2S]- 

As an example, the release of histamine is promoted by orexin, nora- 
drenaline and acetylcholine in the BF and inhibited by GABA in the VLPO. 
Histamine is autoinhibitory, like most of the monoamines. The equation for 
histamine's rate of change, like all of the other neurochemicals involved, also 
requires an adenosine term. Histamine in the human sleep-wake system has 
no known uptake mechanism, and so only degrades as a method of removal. 



13 



Therefore the equation for histamine is 
dH 



dt 



= CiOX+Ci+iNA+Ci+2AChBF+Ci+3AD-Ci+4GABAvLPO-{ci+5-Ci+6)H 



where the q's are rate coefficients. The rest of the equations are formulated 
in the same fashion. 



Fast Time System 

dGABAsFw 



dt 

dGABABFs 



dt 



dOX 
dt 

dH 

~dt 



dAChsF 
dt 

dAChiDT/ppT 



dt 



dNA 
dt 

dS 
~dt 



dDA 
dt 



CiAChBF + C2OX + C3AD - C4GABABFW- 

cr,{GAl){GABABFj (3.2) 
cqAD + C7NA - csGABAbfs - C9{GA1){GABAbfs) 

(3.3) 

CioAChLDT/PPT + (Cll - Ci2)NA + (Ci3 - Ci4)0X + C15AD- 

cieS - cnGABAvLPO + c^^GABAbfs (3.4) 

C19OX + C20NA + C2lAChBF + C22AD - C23GABAVLPO- 
(C24 - C25)H (3.5) 
C260X + C27H + C28AD - C29S + CsoGABAbFs - CsiAChBF 

(3.6) 

C32OX + cssH + C34AD + C35NA - cseS - c^-jAChLDT/ppT 

(3.7) 

casOX + csgAChBF + C40AD - cuGABAvLPO- 
(c42-C43)NA-cu(hNET)(NA) (3.8) 

C45OX + C4qH + a-jNA + a^AD - c^gGABAvLPO- 
{c5o-C5i)S-c,2{hSERT){S) (3.9) 

C53OX + C54S + C^^AChLDT/PPT + C^qNA + 

c^sGABAvLPO - C59 {hDAT) {DA) -cqqDA (3.10) 
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The uptake transporter values are contained in the uptake terms in equa- 
tions (3.2), (3.3), (3.8), (3.9) and (3.10). Their values are found from mi- 
crodialysis experiments with human cerebrospinal fluid (CSF), and are as 
follows t26j: GAl = 1, hNET = 0.457, hSERT = 0.463, hDAT = 1.22. 



3.2.2. Uptake Mechanisms 

The neurotransmitters involved in the sleep-wake system, as with any 
neuronal network, have natural rates of removal from the synaptic gap once 
released. These rates can either be contributed to by degradation, neuro- 
transmitter uptake or a combination of both. All of the neurotransmitters 
involved degrade over time, but only some have neurotransmitter uptake 
mechanisms. We will not be considering any of the smaller scale chemi- 
cal dynamics of the neurotransmitter uptake mechanisms. We care only 
whether a neurotransmitter is removed by an uptake mechanism, and will 
not differentiate between plasma membrane-bound transportation or recep- 
tor uptake. The key take away from the uptake transportation data is that 
only GAB A, noradrenaline, serotonin and dopamine have significant synap- 
tic uptake transport mechanisms in humans. These mechanisms are char- 
acterized by the data given constants GAl, hNET, hSERT and hDAT 
for GABA, noradrenaline, serotonin and dopamine, respectively [26] . These 
neurotransmitter transporters appear in the uptake term of the equations, 
always multiplied by the neurotransmitter it is dependent on to indicate a 
linear removal rate. As an example, hSERT appears in equation 9 multi- 
plied by — 0525*, so as to act as a removal mechanism for serotonin when it is 
present. However, hSERT does not multiply the natural degradation term 

— C^qS. 



3.2.3. Nonlinear Dynamics 

There are two neurotransmitters, AD and GABAvlpo-, operating on a 
slow time scale that define the sleep-wake cycle. The qualitative portrait 
in phase space is known from basic sleep-wake cycle activity. We can con- 
struct the desired closed orbit and concentric curves from the nullclines of 
the AD and GABAvlpo equations. We identify the bounding region in 
the phase space activity by the box contained between {AD, GABAvlpo) = 
{ADmax,GABAmin) sjid {ADmin, GAB Amax) ■ Wc choose small values of 
ADmin = GABAmin = 0.001 rather than exactly zero to prevent inclusion of 
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a potentially unstable steady state. We can scale the size of this box accord- 
ing to the values of ADmax and GABAmin- These two points are represented 
by the intersections of the AD and GABAvlpo nuUclines. To describe the 
activity between these two points, we begin with the fact that when AD is at 
it's maximum value, GABAvlpo is at it's minimum value (0.001). This sig- 
nifies that at the peak of wake activity, sleep activity is at its minimum. As 
we move along the AD nullcline, the concentration of AD decreases to 0.001 
as the concentration of GABAvlpo increases to GABAmax- From the point 
{GABAvlpo, AD) = {GABA„^ax, 0.001), we move along the GABAvlpo 
nullcline slowly until we reach the point {GABAvlpo, ^D) = (0.001, AD^ax) 
again, creating a closed orbit. Figure |8] shows this movement in the phase 
plane, as well as showing the location of the fixed point and neighboring 
trajectories. Everything is in non-dimensional form, and therefore should be 
thought of in terms of threshold levels rather than micro-moles per liter. 

The system is made up of two widely separated time scales: the slow, 24 
hour time scale that is governed by the build-up and decay of adenosine, and 
the fast, 0(1) (minute) time scale that involves interactions of all the sleep 
and wake active neutrotransmitters, but is regulated by the oscillation of the 
slow time scale. 

We begin by forming the slow time system that will later be shown to 
be responsible for the limit cycle dynamics of the system. This system is 
comprised of two equations: the adenosine equation that fluctuates with 
wake activity, and the GABA in the VLPO equation that fluctuates with 
NREM sleep activity. We consider NREM sleep to be the ground state of 
the sleep-wake system, and all other oscillations to be deviations thereof. The 
adenosine equation must involve a constant source term, due to adenosine's 
constant creation as a by-product of neuronal activity, as well as a decay 
term that is dependent on adenosine's natural rate of degradation as well as 
it's loss from the synaptic gap through the its disinhibition of GABA in the 
VLPO [7|. It is well known that mechanism for receptor binding is nonlinear 
[26], and so the disinhibition effect of GABAvlpo via adenosine should be 
nonlinear. It also follows that this nonlinear term should appear in both the 
equations for AD and GABAvlpo because adenosine's disinhibition of the 
VLPO both depletes the concentration of adenosine while simultaneously 
(and with the same rate) adding to the GABAvlpo concentration. The 
GABAvlpo equation must contain the uptake term for GABA, a degrada- 
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(ADmin.GABAmax) 



GABA 



Figure 8: The ideal activity of the two neurotransmitter nuUchnes in the GABAvlpo^AD 
phase plane. 



tion term and a contribution from adenosine as well as the inhibitory and 
excitatory of other neurotransmitters. A constant loss term is required in the 
GABAvLPO equation as a representation of its quick breakdown time in the 
synaptic gap [10]. Though it is being broken down quickly and continuously, 
GABAvLPO still provides a tonic effect on the sleep-wake system as a whole 
|40j . Based on these requirements as well as the observations of the glycolysis 
reaction [33], we obtain the following slow time system: 

Slow Time System 

^ = A;i - k^AD - {GABAvLPO?{AD) (3.11) 



HG A n A 

'"'-'''^ = -e-hGABAvLPo + hAD + iGABAvLPo)\AD) (3.12) 
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where 

e = fi — {aiAChsF + 0'2AChL£iT/ppT + asNA + a^S + a^DA) 

To provide a more physical explanation of the assumed reaction system, we 
know that adenosine is a bi-product of the constant source of energy in the 
nervous system, adenosine triphosphate (ATP). We know that adenosine is 
continuously converted back into ATP |3^. From experiments in the VLPO, 
we know that activation of VLPO neurons through adenosine-mediated inhi- 
bition contributes to the somnegenic actions of adenosine, and only adenosine 
acts directly on the VLPO in the sleep- wake system [3, |3l [301 EH]. During 
wake, increased levels of adenosine promote sleep by disinhibiting the release 
of GABA in VLPO. Adenosine levels in this area increase during wake and 
decrease during sleep 0. In this way the physical foundation for the glycol- 
ysis system is similar to and provides a good framework for our slow time 
scale system. 



3.2.4. Characterizing REM/NREM 

REM and NREM sleep are determined by a nonlinear combination of 
neurotransmitters. Specifically, GABAyipo promotes the state of REM 
sleep, as does acetylcholine from both the LDT/PPT and the BP. While 
NREM can be seen as an absence of the vigilance state and the contributions 
from GABA, REM sleep must incorporate contributions from neurotransmit- 
ters active in the wake state as well as those in the NREM state. Letting 
R(t) represent the REM sleep state, we assume that R" + $-R' + VR = 0. 
To determine <I>, note that the transition from the waking state to the 
sleeping state, and vice versa, occurs when GABAvlpo = AD. Given 
that the REM state has an asymptotically stable steady state during wake, 
which occurs when AD > GABAvlpo^ and has a limit cycle during sleep, 
which occurs when AD < GABAvlpo^ then $ contains a term of the form 
{GABAvlpo - AD){R^ - R^). When entering wake, so GABAvlpo = AD 
and AD > 0, there is an impetus for the impending dormancy of REM and 
the resulting mathematical requirement is that ^{AD, AD) > 0. Likewise, 
when entering sleep, so GABAvlpo = and AD < 0, the requirement for 
the oncoming REM limit cycle is ^{AD, AD) < 0. This gives us 



d'R 



+ 



a{GABAvLPo - AD){R' - R^) + (3 



dAD 
dt 



^ + TR = 0. (3.13) 
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3.2.5. Some Special Cases 



Given the complexity of (3.13), it is worth considering the solution in 



particular situations and then using this to help estimate the oscillatory 
response of equation (3.13). For example, the period for R can be determined 



from the fact that there are three to four normal REM oscillations over a 
normal night's sleep. So, when GABAvlpo = AD, 



(fR dADdR ^„ 

+ 1^— - + TR 



0, 



(3.14) 



df^ ' dt dt 

Given the slowly varying nature of AD., then R = e***, where r = ^{—/3u ± 
fi'^u^ — 4r) and u = (AD)'. Since u strongly controls the decay, and to 
have r control the number of REM oscillations, we require r = —u ± 
where ^/^ is approximately the number of oscillations. Thus, /3 = 2 and 
r = 7 + M^. From this we have that 



d^R 



{GABAvlpo - AD){R'^ - 1.3) + a 



dAD 
dt 



dR 

~di 



+ 7 + 



/dADV 

\dr) 



i? = 0. 



(3.15) 



where a is determined numerically and 8 < 7 < 16. 



Based on Lienard's Theorem, we know that the above oscillator has a 
unique, stable limit cycle [36] • Thus the REM oscillations will continue to 
function in a stable, periodic fashion even when subjected to reasonable per- 
turbations. Evidence of this result is in the Orexin Knockout section. 



3.3. Analysis 

To find the values of the coefficients of the system and to ensure stabil- 
ity, we present several methods of analyzing the equations. We begin with 
the limit cycle coming from the Slow Time System. As stated earlier, the 
adenosine equation drives the wake state dynamics, and the GABA equation 
drives the sleep state dynamics. The specific equations are 

^AD = ki- k2AD - (GABAvlpo? {AD) (3.16) 
dt 

^GABAvlpo = -e-hGABAvLPo+hAD+{GABAvLPo?{AD), (3.17) 
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where 



e = /i — {aiAChsF + ci2AChLi:,T/ppT + OsA^A + a^S + a^DA). 

It is not possible to determine the rate constants with the available ex- 
perimental evidence. So, we take a different tack and establish bounds on 
the constants that will produce the responses known to occur in sleep and 
wake. To begin, we compute the nuUclines of the slow time system in order to 
create a trapping region for the limit cycle and to determine the coefficients 
of the system by using an approximate fixed point. We have four coefficients 
to determine (fci, A;2, fcs, ^4) as well as reasonable values for the parameter 
e. Plugging the approximate fixed point into both nuUclines will give us two 
equations with five unknowns. In order to make the system well determined, 
we observe the maxima and minima of the equations as well as the linearized 
slow time system with the approximate fixed point. Knowing the form of 
the determinant and the trace of this system needed to maintain stability, 
we can determine the rest of the unknowns. 



Now, solving the nuUcline equations obtained from 3.16 and (3.17) for 
AD and we have 

4 P - ^1 

k2 + GABAl^po ' 

and 

^ hGABAvLPo - e 
+ {GABAvLPo? ■ 

Setting these equal to each other, we have three points of intersection, only 
one of which is real. Using GABAvlpo ^ [0,2] and AD G [0,2] as the 
boundaries of the trapping region, we can approximate a fixed point inside 
the closed orbit at GABAvlpo = 0.703 and AD = 0.823. This provides us 
with two of the needed equations. 



We also know that each of these equations is zero separately when AD = 
ADmax and GABAvlpo = GABAmin, as well as when AD = ADmm and 
GABAvlpo = GABAmax- We let AD^ax = GABA^ax = 2 and AD^in = 
GABAmin = 0.01. Taking e = 0.3, and using the stated values for the 
max and min concentrations it follows that ki = 0.49, /c2 = 0.1, k^ = 0.3. 
^4 = 0.15. We also note that 0.29 < e < 0.32 for the system to remain stable. 
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3.3.1. Remaining Coefficient Estimation 

Now that the coefficients for the slow time system are determined, we can 
use the requirements that the eigenvalues have either pure imaginary form 
or negative real parts in the complex form to determine the coefficients of 
the entire system. We have the solution vector of the linearized system as x = 
[GABAbfw, GABAbfs, ox, H, AChBF, ACHldt/ppt, NA, S, DA, AD, GAB A 
The matrix of coefficients and linearized values for the system are given in 
Table El 

In order for the system to be stable, we require that the real parts of the 
eigenvalues to be negative. Restricting the coefficients to all be between 
and 1, we use an algorithm that generates coefficients that keep the real parts 
of the eigenvalues negative. We also note that the coefficients in GABAvlpo 
equation need to be small compared to the rest of the coefficients in order to 
keep e bounded. 

3.3.2. Sub-Systems 

When analyzing such a large system of equations with undetermined co- 
efficients, it is helpful to break the large system into smaller sub-systems and 
observe the dynamics on a local scale. We can utilize the neurophysiology 
of the physical neurotransmitter sub-systems to determine the coefficients as 
well as the stability proffies of the individual sub-systems. We break the sys- 
tem up into to following categories: NREM Specific, Sleep Specific (NREM 
Inclusive), Wake Specific, REM + NREM Specific, REM + Wake Specific, 
Wake Inclusive, Sleep Inclusive, REM Inclusive. Note that there is no set 
of equations associated with REM only. The categories represent what neu- 
ron populations are active during the corresponding states. We observe the 
steady-state behavior of each sub-system. Since we know the relative mag- 
nitude of the neurotransmitter amounts as well as their dynamic behavior 
during the states of REM, NREM and Wake, the coefficients of these smaller, 
more manageable systems should be more easily attainable. This also allows 
us to observe the more manageable systems for potentially extreme behavior, 
and provides us with bounds to prevent this from happening. This analy- 
sis does not give us specific values for the system coefficients, knowing the 
relative values makes visualizing the system easier as well as simplifying the 
eigenvalue analysis. 
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NREM Specific 

The NREM specific system involves the neuron population only active 
during the state of NREM sleep. The neuron population that is only active 
in NREM sleep is the GABAergic population in the BP. We assume all 
equations but ''■'^^^^^bfs |-,g equal to zero, making each neurotransmitter 
other than GABAbfs constant, and assume the other neurotransmitters in 
the equation to be at their steady state values when only GABAbfs is active. 
This is not necessarily a physiologically accurate view of the system for an 
extended length of time, but we can view a snapshot of this activity. We 
follow this same criteria for the rest of the small systems. 

Wake Specific 

The Wake specific system involves neuron populations with the most 
activity during the wake state. The monoaminergic neuron populations (in- 
cluding OX, H, NA, S and DA) are most active in this state. 

Sleep Specific (NREM Inclusive) 

The Sleep Specific system involves neurotransmitters that promote sleep 
states (REMS and NREMS). The populations that are active in this state 
are GABAergic populations in the VLPO and BF. 

REM + Wake Specific 

The REM and Wake Specific system involves neurotransmitters that pro- 
mote REMS and the Wake state. The neuron populations that are active 
during these states are GABAergic populations in the BF and chohnergic 
populations in the BF and LDT/PPT. 

REM + NREM Specific 

The REM and NREM Specific system involves neurotransmitters promot- 
ing REM sleep and the NREM sleep. The neuron population that is active 
during these states is the GABAergic population in the VLPO. 

Wake Inclusive 

The Wake Inclusive system involves neurotransmitters that promote the 
Wake state as well as those that promote both Wake and REMS. The neuron 
populations that are active during these states are the GABAergic population 
in the BF, the cholinergic population in the BF and LDT/PPT, and the 
monoaminergic populations. 
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Sleep Inclusive 

The Sleep inclusive system involves neurotransmitters associated with 
NREM specific, REM and NREM specific, and Wake and REM specific 
states. The neuron populations that are active during these states are 
GABAergic populations in the VLPO and BP, and cholinergic populations 
in the BP and LDT/PPT. 

REM Inclusive 

The REM Inclusive system involves neurotransmitters associated with the 
Wake and REM specific and REM and NREM specific states. The neuron 
populations that are active during these states are the GABAergic popula- 
tion in the VLPO, the wake promoting GABAergic population in the BP as 
well as the cholinergic populations in the BP and LDT/PPT. 



This analysis allows for an approximation of the coefficients, but it does 
not yield exact values for these coefficients. Pigures 11 and 12 show the 
change in concentration levels of Acetylcholine in the BP and GABA in the 
BP, respectively over a period of 24 hours. The remainder of the neurotrans- 
mitters involved show the same fiuctuations as that of Acetylcholine, whereas 
GABA in the sleep active portion of the BP is unique. 




Figure 9: Concentration levels of 
Acetylcholine in the BF over 35 



CihouJ) " " " hours. 



4. Results and Discussion 

The results of the implemented normal, unperturbed, system are shown 
in Pigures 11 14 Pigure 11 depicts the slow time variables (which are the 
easiest to visualize on a large time scale) AD and GABAvlpo- The AD 
curve represents the non-dimensional amount of adenosine accumulation in 
the BP, which consequently measures sleep pressure. The GABA curve rep- 
resents the non-dimensional amount of GABA accumulation in the VLPO, 
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Figure 10: Concentration levels 
of GABA in the BF associated 
with sleep over 35 hours. 



which is a measure of wake pressure. Both exhibit a rapid transition be- 
tween sleep and wake states as well as ongoing periodic behavior for nine 
days. The system runs like this for much larger periods, demostrating that 
unless perturbed, the sleep-wake cycle is a self-sustained, periodic oscillation. 
The period of each cycle is approximately 24 hours - the wake state occupies 
roughly two-thirds of this time, while the sleep state makes up the remain- 
ing one-third. This is consistent with our knowledge of the structure of the 
human sleep-wake cycle. Figure 12 depicts one cycle in the GABA — AD 
phase plane, showing the exact behavior we required from our previous nuU- 



cline analysis. Figure 13 is the isolated REM curve for 30 hours. Figure 



14 shows the AD, GABA and REM curves simultaneously in a sample 48 



hour time slice. The REM/NREM oscillations respond quickly to the tran- 
sitions between wake and sleep, and cycle through three to four oscillations 
in the designated eight hour sleep period. From the biological background 
and neurotransmitter details, this displays physiologically reasonable results. 



4.1. Orexin Knockout Experiment 

Narcolepsy is a neurodegenerative disease where there is a loss of neu- 
rons containing orexin (OX) [13]. It has been shown to exist across many 
mammals, and has been studied extensively in mice. Here, we take infor- 
mation from the sleep-wake cycles of orexin knock-out (which means the 
targeted destruction of orexin neurons) mice. These mice exhibit the nar- 
colepsy symptoms of sleep fragmentation, sleep onset during REM sleep 
states and increased spontaneous REM and NREM sleep during the vigi- 
lance state [271 SD H ESI Eg. 

A reasonable method to represent the orexin knock-out perturbation to 
the system is to change the flux of the OX amount equation. Rather than 
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Figure 11: Adenosine and GABA (in the VLPO) levels of the sleep-wake cycle model over 
216 hours (9 days). 



take orexin out all together, we limit the rate at which orexin is produced, 
but keep the removal factors the same. This is a reasonable way to represent 
the physiological destruction of orexin neurons, since none of the uptake and 
removal mechanisms are necessarily effected. It has been shown that levels of 
orexin decrease in human cerebral spinal fluid by at least 80 percent |13], and 
so we multiply orexin producing elements of the OX equation {AChLDT/ppT, 
NA, OX and AD) by 0.2. Implementing these changes numerically, the 
periodicity and limit cycle of the system is maintained with orexin decrease. 



This is illustrated in Figure 15, which clearly shows the same characteristics 
as seen in Figure 11 However, we also don't see the stereotypical changes 



to the transition rates between behavioral states characteristic to narcolepsy. 
We have shown that the system can withstand perturbations [27], but may 
need to think of implementing another method of making genetically specific 
perturbations. This also shows that orexin may have to play a larger role 
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Figure 12: Phase plane of sleep-wake cycle model for one cycle. 



in measuring behavioral state instability, and may have to withstand more 
non-linear changes to its equation when investigated further. 

4-2. Sleep Camp Data 

[6] describes experiments related to perturbations in an individuals sleep/wake 
cycle. The trials performed were carried out in a controlled situation with 
a group of 40 adolescents over a period of one week. The perturbations 
performed were delayed sleep and wake onset, prolonged sleep and wake pe- 
riods and forced wake. These perturbations were achieved without the use of 
external chemical stimuli, and so were essentially achieved in a noninvasive 
fashion. Salivary Cortisol and melatonin levels were recorded daily as well 
as frequent EEG readings for accurate brainwave activity. The initial time 
and length of the perturbation period were recorded, and so this data was 
easy to translate into perturbations in the model. About 50 perturbations 
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Figure 13: REM cycle visualization of sleep-wake cycle model over a day. 



were introduced (which translates to a period of about 10 days) and their 
effect on the system was recorded using graphical and data analysis methods. 



The perturbations to the system were only introduced to one variable 
- adenosine. Forced sleep and wake were achieved by forcing AD to its 
maximum at the appropriate forced waking initialization, and forcing AD to 
its minimum at the appropriate forced sleep time. These sleep and wake times 
as well as the prolonged sleep and wake periods were taken from a data set of 
seven adolescent individuals. Once perturbed the amount of times done so in 
the physical experiment, we allow the system to run without perturbation for 
several cycles. This entire set of perturbations and the response is depicted in 



Figure [T6| Figure [T7| shows what the subject's cycle would look like without 
perturbation. One can see that the perturbations cause small scale (on the 
order of a day or two) changes in the length of sleep or wake state. However, 
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Figure 14: Unperturbed system of AD, GABA and REM from hours 150 to 200 - over two 
days. 



once allowed to run unperturbed, the system settles back into a regular 24 
hour period. Initial wake and sleep times have drifted, but that is to be 
expected. When run for a longer time, the system will eventually drift back 
into its original phase with the same wake and sleep initialization times as the 
subject started with. We know this is an accurate response to these non-drug 



related, small perturbations(REF). Figures 18 and 19 show a quantitative 
measure of the changes in period in the unperturbed vs. the perturbed 
cycles. Figure 20 demonstrates the upholding of limit cycle dynamics in the 



face of perturbations, and also that the AD — GABA phase plane retains its 
size and shape but drifts in phase. 



5. Conclusions 

We have presented a model of the human sleep-wake generating system 
that identifies specific neuronal regions in the brain with their associated neu- 
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Figure 15: Orexin knockout system over 216 hours (9 days). 



retransmitters and projections. It accounts for the homeostatic influences of 
sleep-wake, and has the potential to encompass circadian influences as well. 
It improves upon previous models of sleep-wake cycles in that the oscilla- 
tions are self-sustaining without the use of external periodic forcing, and it 
is formulated specially to look at the specific sleep-wake cycle of humans. 
We use the known neurochemical dynamics of the human sleep-wake cycle to 
derive equations and form the stability criteria for mathematical phenomena. 
Since based on limit cycle dynamics, the system runs with the contributions 
of all the necessary neurotransmitters and with no need of an external driv- 
ing force. This system gives reasonable predictions for "natural" sleep bout 
duration, wake bout duration, REM bout duration and NREM bout dura- 
tion, as well as estimated natural sleep initiation and wake initiation times. 
The NREM sleep state is considered the ground state from which all other 
activity deviates, and REM activity is determined from the difference of the 
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Figure 16: Perturbed slow time system with forced wake onset and forced sleep onset over 
a one week period. Extended time to show phase drift and maintained periodicity. 



wake and NREM states, as well as the contributions from REM active neuro- 
transmitters. Using experimental data from a mammalian orexin knockout 
system as well as simply and non-chemically perturbed human sleep-wake 
cycles, we were able to visualize the perturbed cycle and compare to the 
experimental results without instability occurring. We can draw a number 
of conclusions from the resulting model. In the human systems, given forced 
prolonged waking, sleep duration reacts accordingly to build up appropriate 
wake pressure. Once the cycle is restored to normal sleep onset and wake on- 
set, the same period for an unperturbed system is maintained. The system 
exhibits some phase drift with this forcing. Given forced prolonged sleep, 
wake duration also reacts accordingly to build up appropriate sleep pressure. 
Once the cycle is restored to normal sleep and wake onset, the periodicity is 
maintained. The system exhibits some phase drift with this forcing. Given 
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Figure 17: Unperturbed cycle allowed to run for the same amount of time as the perturbed 
cycle. 



forced wake up during sleep, the subject does not stay awake for long, and 
falls back into a normal sleep schedule, yielding only some phase drift that 
is corrected after a few cycles. A subject's REM/NREM cycle goes through 
three or four oscillations per state (REM or NREM). The first NREM sleep 
stage is triggered by sleep onset, or when adenosine levels fall below a cer- 
tain threshold and GABA becomes the dominant neurotransmitter. The last 
REM cycle is turned off by wake onset, or when adenosine pressure builds 
back up and GABA decays. Subjects naturally fall asleep and wake up at the 
intersection points of the GABA and Adenosine thresholds - falling asleep 
when enough sleep pressure from AD has been accumulated, and waking up 
when the sleep pressure has decayed once again and when the subject has 
gone though three to four cycles each of REM and NREM sleep. We can see 
that given no external forcing (no stimuli, drugs, etc), a "normal" human's 
sleep- wake schedule unfolds in a periodic, 24 hour fashion that endures in- 
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Figure 18: Examination of the length one period in an unperturbed cycle. 



definitely unless given some sort of perturbation. 

It is important to note the implications of the differences between home- 
ostatic and circadian drive, and the importance of both of these drives' effect 
on sleep. Unlike most recent neurochemically driven models, ours has no 
need for external forcing to keep it periodically driven. The periodic, limit 
cycle activity arises naturally from the neurotransmitter dynamics and the 
equations themselves. This demonstrates an intrinsic, homeostatic force that 
keeps the human sleep-wake cycle functioning without a change in periodicity, 
regardless of whether circadian control is included or not. The perturbations 
in the experimental section demonstrate that even when subjected to strong 
forcing, the sleep-wake cycle retains its natural 24 hour period, and shifts in 
phase are the only real artifacts of outside forcing. The next step would be 
to develop a more complete predictive model of the human sleep wake cycle 
for use in the scientific and medical fields. 
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Figure 19: Examination of the length of one period in a perturbed cycle. 
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Wake 


Inactive 


Inactive 


Active 


High firing 
rate, high 
level. 


Tonic, low 
firing rate, 
high level. 


High firing 
rate, high 
level. 


Fast 

cortical 

rhythms 


Tonic, low 
firing rate, 
high level. 


Tonic, low 
firing rate, 
high level. 


No rate 
change. 
Elevated 
levels. 


Levels in- 
crease with 
prolonged 
wake in the 
BF. 


Projects to: 


LH, LC, 
LDT/PPT, 
vPAG, RN, BF 


BF, VLPO 


BF, VLPO 


LH, TMN, RN, 
BF 


TMN, LH 


BF, VLPO 


VLPO 


LC, BF, RN, 
VLPO 


RN, VLPO, LH, 
BF 


vPAG, VLPO 


Not Applicable 


Receives Projec- 
tions From: 


LC, vPAG, RN, 
BF 


BF, RN, LC, LH 


BF, RN, LC, LH 


LH, VLPO, RN, 
TMN 


TMN, LH, 
VLPO 


nBF, RN, LC, LH 


nBF, VLPO 


O 

> 


RN, LH, LC, 
VLPO 


vPAG, VLPO 


Not Applicable 


1 Neurotrans 
-mitter 


) GABA 


GABA 


GABA 


Orexin 


Histamine 


Acetylcholi 


Acetylchoh 


Noradrenal 


Serotonin 


Dopamine 


Adenosine 
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Parameter 


Value 


Parameter 


Value 


GAl 


1 


C27 


0.10973 


hNET 


0.457 


C28 


0.32943 


hSERT 


0.463 


C29 


0.57879 


hDAT 


1.22 


cao 


1 


ki 


0.49 


C31 


0.02091 


k2 


0.1 


C32 


0.12648 


h 


0.3 


C33 


0.23472 




0.15 


C34 


0.57122 




0.3 


C35 


0.02332 


a 


1 


C36 


1 


7 


8 


C37 


0.61305 


Cl 


0.75709 


C38 


0.06864 


C2 


0.28014 


C39 


0.08638 


C3 


0.61048 


C40 


A 1 

0.1 


C4 


0.76636 


C41 


0.543 


C5 


0.32431 


C42 





C6 


0.83153 


C43 


0.13822 


cy 


0.03471 


C44 


0.35956 


C8 


0.01 


C45 


0.11839 


cg 


0.19577 


C46 


0.13753 


ClO 


0.79157 


C47 


0.1 


Cll 


0.97026 


C48 


0.537 


Cl2 


-1 


C49 





Cl3 





C50 


0.67749 


Ci4 


0.36341 


C51 


0.53609 


Cl5 


0.70633 


C52 


0.36464 


Cl6 


0.1 


C53 


0.59591 


Cl7 


0.97643 


C54 


0.75091 


Cl8 


0.56740 


C55 


A 1 

0.1 


Cl9 


0.91859 


C56 


0.22 


C20 


0.50364 


C57 


10-5 


C21 


0.23758 


C58 


10-5 


C22 


0.1 


C59 


10-5 


C23 


1 


C60 


10-5 


C24 





C61 


10-5 


C25 


0.92037 






C26 


0.04185 







Table 3: Coefficient and parameter values used in the sleep-wake cycle model. 
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Figure 20: Perturbed GABA- Adenosine phase plane. 
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